



multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) {
  library(grid)

  # Make a list from the ... arguments and plotlist
  plots <- c(list(...), plotlist)

  numPlots = length(plots)

  # If layout is NULL, then use 'cols' to determine layout
  if (is.null(layout)) {
    # Make the panel
    # ncol: Number of columns of plots
    # nrow: Number of rows needed, calculated from # of cols
    layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
                    ncol = cols, nrow = ceiling(numPlots/cols))
  }

 if (numPlots==1) {
    print(plots[[1]])

  } else {
    # Set up the page
    grid.newpage()
    pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))

    # Make each plot, in the correct location
    for (i in 1:numPlots) {
      # Get the i,j matrix positions of the regions that contain this subplot
      matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))

      print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
                                      layout.pos.col = matchidx$col))
    }
  }
}    

library(interplot)

model1 <- glm(purge3 ~ dp1 + medianp1 + dp1*defense + defense + logmilex  +  firstleader + IS_WAR  +lrgdpch + mildic + royal + party + pastpurges3+ decade1970s + decade1980s + decade1990s + decade2000s + purge3yrs +
                   purge3yrs2 + purge3yrs3  ,family = binomial (link="logit"))

defplot = interplot(m = model1, var1 = "defense", var2 = "dp1", hist=TRUE) + 
    # Add labels for X and Y axes
    xlab("") +
    ylab("Estimated marginal effect on Pr(purge)") +
    # Change the background
    theme_bw() +
    # Add the title
    ggtitle("Defense alliances") +
    theme(plot.title = element_text(size = 20, face = "bold")) +
    theme(plot.title = element_text(face="bold")) +
    theme(axis.text=element_text(size=12),
        axis.title=element_text(size=10,face="bold")) +
    # Add a horizontal line at y = 0
    geom_hline(yintercept = 0, linetype = "dashed")        
    
    
model2 <- glm(purge3 ~ dp1 + medianp1 + dp1*consultneutralnonaggonly + consultneutralnonaggonly + logmilex  +  firstleader + IS_WAR  +lrgdpch + mildic + royal + party + pastpurges3 + decade1970s + decade1980s + decade1990s + decade2000s + purge3yrs +
                   purge3yrs2 + purge3yrs3  ,family = binomial (link="logit")) 
                   
neutralplot = interplot(m = model2, var1 = "consultneutralnonaggonly", var2 = "dp1", hist=TRUE) + 
    # Add labels for X and Y axes
    xlab("") +
    ylab("") +
    # Change the background
    theme_bw() +
    # Add the title
    ggtitle("Neutrality agreements") +
    theme(plot.title = element_text(size = 20, face = "bold")) +
    theme(plot.title = element_text(face="bold")) +
    theme(axis.text=element_text(size=12),
        axis.title=element_text(size=10,face="bold")) +
    # Add a horizontal line at y = 0
    geom_hline(yintercept = 0, linetype = "dashed")            
    
model3 <- glm(purge3 ~ dp1 + medianp1 + dp1*lsecdum + lsecdum + logmilex  +  firstleader + IS_WAR  +lrgdpch + mildic + royal + party + pastpurges3 + decade1970s + decade1980s + decade1990s + decade2000s + purge3yrs +
                   purge3yrs2 + purge3yrs3  ,family = binomial (link="logit")) 
                   
secdumplot = interplot(m = model3, var1 = "lsecdum", var2 = "dp1", hist=TRUE) + 
    # Add labels for X and Y axes
    xlab("Elite strength") +
    ylab("Estimated marginal effect on Pr(purge)") +
    # Change the background
    theme_bw() +
    # Add the title
    ggtitle("U.S. security assistance") +
    theme(plot.title = element_text(size = 20, face = "bold")) +
    theme(plot.title = element_text(face="bold")) +
    theme(axis.text=element_text(size=12),
        axis.title=element_text(size=10,face="bold")) +
    # Add a horizontal line at y = 0
    geom_hline(yintercept = 0, linetype = "dashed")              
    
model4 <- glm(purge3 ~ dp1 + medianp1 + dp1*larmsdum + larmsdum + logmilex  +  firstleader + IS_WAR  +lrgdpch + mildic + royal + party   + pastpurges3+ decade1970s + decade1980s + decade1990s + decade2000s + purge3yrs +
                   purge3yrs2 + purge3yrs3  ,family = binomial (link="logit")) 
                   
armsplot = interplot(m = model4, var1 = "larmsdum", var2 = "dp1", hist=TRUE) + 
    # Add labels for X and Y axes
    xlab("Elite strength") +
    ylab("") +
    # Change the background
    theme_bw() +
    # Add the title
    ggtitle("Arms shipments (SIPRI)") +
    theme(plot.title = element_text(size = 20, face = "bold")) +
    theme(plot.title = element_text(face="bold")) +
    theme(axis.text=element_text(size=12),
        axis.title=element_text(size=10,face="bold")) +
    # Add a horizontal line at y = 0
    geom_hline(yintercept = 0, linetype = "dashed")        
    
      
multiplot(defplot, secdumplot, neutralplot, armsplot, cols=2)      




    
    

